home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / DMFSD.FOR < prev    next >
Text File  |  1985-11-29  |  4KB  |  112 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DMFSD
  5. C
  6. C        PURPOSE
  7. C           FACTOR A GIVEN SYMMETRIC POSITIVE DEFINITE MATRIX
  8. C
  9. C        USAGE
  10. C           CALL DMFSD(A,N,EPS,IER)
  11. C
  12. C        DESCRIPTION OF PARAMETERS
  13. C           A      - DOUBLE PRECISION UPPER TRIANGULAR PART OF GIVEN
  14. C                    SYMMETRIC POSITIVE DEFINITE N BY N COEFFICIENT
  15. C                    MATRIX.
  16. C                    ON RETURN A CONTAINS THE RESULTANT UPPER
  17. C                    TRIANGULAR MATRIX IN DOUBLE PRECISION.
  18. C           N      - THE NUMBER OF ROWS (COLUMNS) IN GIVEN MATRIX.
  19. C           EPS    - SINGLE PRECISION INPUT CONSTANT WHICH IS USED
  20. C                    AS RELATIVE TOLERANCE FOR TEST ON LOSS OF
  21. C                    SIGNIFICANCE.
  22. C           IER    - RESULTING ERROR PARAMETER CODED AS FOLLOWS
  23. C                    IER=0  - NO ERROR
  24. C                    IER=-1 - NO RESULT BECAUSE OF WRONG INPUT PARAME-
  25. C                             TER N OR BECAUSE SOME RADICAND IS NON-
  26. C                             POSITIVE (MATRIX A IS NOT POSITIVE
  27. C                             DEFINITE, POSSIBLY DUE TO LOSS OF SIGNI-
  28. C                             FICANCE)
  29. C                    IER=K  - WARNING WHICH INDICATES LOSS OF SIGNIFI-
  30. C                             CANCE. THE RADICAND FORMED AT FACTORIZA-
  31. C                             TION STEP K+1 WAS STILL POSITIVE BUT NO
  32. C                             LONGER GREATER THAN ABS(EPS*A(K+1,K+1)).
  33. C
  34. C        REMARKS
  35. C           THE UPPER TRIANGULAR PART OF GIVEN MATRIX IS ASSUMED TO BE
  36. C           STORED COLUMNWISE IN N*(N+1)/2 SUCCESSIVE STORAGE LOCATIONS.
  37. C           IN THE SAME STORAGE LOCATIONS THE RESULTING UPPER TRIANGU-
  38. C           LAR MATRIX IS STORED COLUMNWISE TOO.
  39. C           THE PROCEDURE GIVES RESULTS IF N IS GREATER THAN 0 AND ALL
  40. C           CALCULATED RADICANDS ARE POSITIVE.
  41. C           THE PRODUCT OF RETURNED DIAGONAL TERMS IS EQUAL TO THE
  42. C           SQUARE-ROOT OF THE DETERMINANT OF THE GIVEN MATRIX.
  43. C
  44. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  45. C           NONE
  46. C
  47. C        METHOD
  48. C           SOLUTION IS DONE USING THE SQUARE-ROOT METHOD OF CHOLESKY.
  49. C           THE GIVEN MATRIX IS REPRESENTED AS PRODUCT OF TWO TRIANGULAR
  50. C           MATRICES, WHERE THE LEFT HAND FACTOR IS THE TRANSPOSE OF
  51. C           THE RETURNED RIGHT HAND FACTOR.
  52. C
  53. C     ..................................................................
  54. C
  55.       SUBROUTINE DMFSD(A,N,EPS,IER)
  56. C
  57. C
  58.       DIMENSION A(1)
  59.       DOUBLE PRECISION DPIV,DSUM,A
  60. C
  61. C        TEST ON WRONG INPUT PARAMETER N
  62.       IF(N-1) 12,1,1
  63.     1 IER=0
  64. C
  65. C        INITIALIZE DIAGONAL-LOOP
  66.       KPIV=0
  67.       DO 11 K=1,N
  68.       KPIV=KPIV+K
  69.       IND=KPIV
  70.       LEND=K-1
  71. C
  72. C        CALCULATE TOLERANCE
  73.       TOL=ABS(EPS*SNGL(A(KPIV)))
  74. C
  75. C        START FACTORIZATION-LOOP OVER K-TH ROW
  76.       DO 11 I=K,N
  77.       DSUM=0.D0
  78.       IF(LEND) 2,4,2
  79. C
  80. C        START INNER LOOP
  81.     2 DO 3 L=1,LEND
  82.       LANF=KPIV-L
  83.       LIND=IND-L
  84.     3 DSUM=DSUM+A(LANF)*A(LIND)
  85. C        END OF INNER LOOP
  86. C
  87. C        TRANSFORM ELEMENT A(IND)
  88.     4 DSUM=A(IND)-DSUM
  89.       IF(I-K) 10,5,10
  90. C
  91. C        TEST FOR NEGATIVE PIVOT ELEMENT AND FOR LOSS OF SIGNIFICANCE
  92.     5 IF(SNGL(DSUM)-TOL) 6,6,9
  93.     6 IF(DSUM) 12,12,7
  94.     7 IF(IER) 8,8,9
  95.     8 IER=K-1
  96. C
  97. C        COMPUTE PIVOT ELEMENT
  98.     9 DPIV=DSQRT(DSUM)
  99.       A(KPIV)=DPIV
  100.       DPIV=1.D0/DPIV
  101.       GO TO 11
  102. C
  103. C        CALCULATE TERMS IN ROW
  104.    10 A(IND)=DSUM*DPIV
  105.    11 IND=IND+I
  106. C        END OF DIAGONAL-LOOP
  107. C
  108.       RETURN
  109.    12 IER=-1
  110.       RETURN
  111.       END
  112.